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Abstract The Landauer approach provides a conceptually simple way to calculate the 

intrinsic thermoelectric (TE) parameters of materials from the ballistic to the diffusive 
transport regime. This method relies on the calculation of the number of propagating 
modes and the scattering rate for each mode. The modes are calculated from the energy 
dispersion (E(k)) of the materials which require heavy computation and often supply 
energy relation on sparse momentum (k) grids. Here an efficient method to calculate 
the distribution of modes (DOM) from a given E(k) relationship is presented. The 
main features of this algorithm are, (i) its ability to work on sparse dispersion data, 
and (ii) creation of an energy grid for the DOM that is almost independent of the 
dispersion data therefore allowing for efficient and fast calculation of TE parameters. 
The effect of k-grid sparsity on the compute time for DOM and on the sensitivity of the 
calculated TE results are provided. The algorithm calculates the TE parameters within 
5% accuracy when the K-grid sparsity is increased up to 60% for all the dimensions 
(3D, 2D and fD). The time taken for the DOM calculation is strongly influenced 
by the transverse K density (K perpendicular to transport direction) but is almost 
independent of the transport K density (along the transport direction). The DOM and 
TE results from the algorithm are bench-marked with, (i) analytical calculations for 
parabolic bands, and (ii) realistic electronic and phonon results for Bi2Tes- 

Keywords Landauers method • Thermoelectricity • Electronic structure • Phonons • 
Density of Modes 



1 Introduction 

Thermoelectricity provides an attractive and a clean way of converting waste heat into 
electricity. There have been a lot of efforts to improve the efficiency of thermoelectric 
(TE) devices. Solid-state TE devices are aggressively pursued both in the industry 
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and research due to their advantages such as, (i) compactness, (ii) resistance to wear 
and tear, and (iii) portabihty. Thermoelectric efficiency (ZT) improvements need very 
careful engineering designs and optimization in terms of, (i) materials [T][51|31|l] , (ii) 
structures like superlattices, nanocomposities, etc. [5l|6l[7l[8] , and (iii) devices [Ol llOlfTTl 
I12lll3j . With so many design parameters it is extremely difficult to experimentally test 
every possible combination. At this point computer modeling plays a very significant 
role in designing and optimizing TE devices from material to the system level [14II15I 
I16| . The present work focuses on the calculation of the TE transport parameters using 
the material energy dispersion as shown in Fig. [T] 

In the present work we focus on the calculation of the material properties involved 
in the calculation of ZT. The value of ZT and the thermoelectric power-factor (PF) 
for a material are given by |17lll8lfT9] . 



where G, S, He and k; are the electronic conductivity, electronic Seebeck coefficient, 
electronic thermal conductivity and lattice thermal conductivity, respectively. The term 
'd' is the dimensionality of the conductor. All the TE parameters depend on the elec- 
tronic and lattice properties of the material. These material properties are strongly 
coupled and an improvement in one of the coefficients may degrade the other [19) . 

The Boltzmann transport equation (BTE) [191120] has been the most commonly 
used method to calculate the TE material parameters. However, with reduced dimen- 
sionality of the TE materials (like nanodots, nanowires, etc) the application of the 
Landauer approach [211122] for calculating the TE parameters has gained a lot of at- 
tention [17lll8l[23ll24ll25ll26ll27] due to the simplicity of the approach. The Landauer 
approach is applicable from the ballistic to the diffusive regime of transport for nano- 
structures. This model is insightful for understanding the impact of dimensionality on 
TE parameters too [18] . 

At the core of the Landauer approach is the calculation of distribution of modes 
(DOM) (17 ^.23^18 which is similar to the transport distribution function (TDF) used 
in the BTE model [28] as shown in Refs. [171123] for both electrons and phonons. The 
DOM represents the number of conducting channels available for the carriers, like elec- 
trons or phonons, at a given energy. From computational aspect, most of the previous 
work using Landauer's approach relied heavily on very fine E(k) calculations and then 
calculating the DOM by band-counting method [l^ as shown in Fig. [2] The BTE 
methods use the reduced Brillouin Zone (BZ) integration schemes [291130] to calculate 
the TE parameters. However, these approaches too depend on a fine momentum mesh 
for numerical integrations. The computation of the dispersion relations in novel ma- 
terials require significant computational resources and in general delivers results on 
momentum meshes that are not dense enough to derive a complete DOM or TDF. 

To overcome the above mentioned computational challenges an efficient algorithm 
to calculate the DOM (used in the Landauer model) from a given E(k) is outlined in 
this work. The present method has two advantages over the previous band-counting 
methods, which are, (i) the energy dispersion (E(k)) can be relatively sparse, and (ii) 
the energy grid for the DOM and the E(k) does not have to be identical. Overall com- 
pute time for the calculation of TE parameters is reduced in two steps, (i) relatively 




(1) 



PF = G-S'^ [W/K'^rn-~'^] 



(2) 
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little compute time is needed to calculate the DOM from the sparse energy disper- 
sion, and (ii) the sparse DOM energy grid further reduces calculation time of the TE 
parameters (G, S, Ke, and K/). 

The present work is divided in the following sections. The basic TE theory in 
the linear transport regime is outlined in Sec. |2.1[ The generic algorithm for DOM 
calculation is presented in Sec. 2.2| with specific changes required for electrons in Sec. 
|2.3[ and for phonons in Sec. |2.4 Discussion on the transmission calculation is provided 
in Sec. |2.5[ The results section provides the k dependent sensitivity analysis in Sec. |3.1| 
and the timing analysis of the algorithm in Sec. |3.2| Comparison and verification of the 
TE parameters, calculated using the algorithm, with published results are provided in 
Sec. |3.4[ The summary of the work is outlined in Sec. |4] 



2 Theory and methodology 

In this section the calculation of the TE parameters and the details of the algorithms 
are outlined. 



2.1 Thermoelectric parameters in the linear transport regime 

The ZT of a material at a temperature, T, is based on the calculation of the intrinsic 
material properties which include both the electronic and the lattice properties (see Eq. 
([T])). The electronic transport parameters are obtained using the Landauer approach 
in the zero current limit ;17..23) as, 

G=^-/o [n-^Tu"-^] (3) 

S = -[kB/q]-[h/Io] [V/K] (4) 
Tk 
li 



^.-j2^;[^lnE).MiE).'-^.,E. (6) 

where Ij is the jth order energy moment integration around the Fermi Level (Ep). 
The terms q, ks and h are the electronic charge, Boltzmann constant, and Planck's 
constant, respectively. In the quantity Ij the terms A4(E), T(E), and J-po are the 
distribution of modes (DOM) at energy E, transmission at energy E and the Fermi- 
Dirac distribution function, respectively. 

The lattice thermal conductivity (k;) can be calculated using the Landauer's model 
as PPSPT], 



^i{T)^hPi [W/m'^-^K] (7) 
Pj = / T{u}) ■ M{uj) ■ J ■ (8) 
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Table 1 Dimensionality of structure and dependence on 'K' vectors 



Structure 


Periodic 


Confined 




i^ll 


K 




(dimension) 


dim. (P) 


dim. (C) 










Wires (ID) 


1 


2 





1 


[^11 




Films (2D) 


2 


1 


1 


1 






Bulk (3D) 


3 





2 


1 







where Pj is the jth order phonon energy integration. The terms M{uj) and T{y>) 
are the distribution of modes, and the transmission of the modes at a phonon frequency 
ui, respectively. 

Equations Q and ([7| show that the calculation of any transport parameter within 
the Landauer model depend on two quantities, (i) the distribution of modes {M) 
and, (ii) the transmission of the modes (T). The DOM depends only on the energy 
dispersion of the carriers in the material whereas the transmission is controlled by the 
dispersion and the scattering mechanisms of the carrier. The advantage of the Landauer 
model lies in the separation of the transport kernel into two parts that can be solved 
using parallel computer programming leading to a faster and efficient calculation of 
the transport parameters. In the next part the details of the algorithm to efficiently 
calculate the DOM from a given energy dispersion is outlined. 



2.2 Distribution of Modes calculation 

The step by step procedure for the calculation of DOM (applicable to both electrons 
and phonons) is given below, 

1. Obtain the energy dispersion of an d dimensional, where d = 1, 2 or 3, periodic 
material. The momentum vector ' K' can be decomposed into two components, (i) 
along the transport direction denoted by A'n , and (ii) in the direction perpendicular 
to the transport direction denoted by K^^ depending on the dimensionality of the 
conductor as shown in Table [T] 

2. For each combination of Kj^, a ID _E — i^y is obtained which is used for mode 
counting. The energy grid for the DOM (EGD) is created based on the ID E — K]^^ 
for all the K± . This energy grid does not have to be identical to the energy values 
from the E{k) data. The details of choosing the energy limits for the electrons and 
phonons are outlined in Sec |2.3| and Sec. |2.4| respectively. The energy grid is chosen 
so as to provide a reasonable compromise between the computation time and the 
accuracy of the results. 

3. For a ID E — , the group velocity (vgrp) is calculated to find out the regions of 
monotonic variation in the energy with K^^ . Only positive A"|| are considered since 
the E — Ku relations are symmetric. The +ve half group velocity is calculated as. 



1 dE{K) 
h dKu 



"arp - (9) 



The monotonic velocity regions, Rl and R2, for an example E(k) are shown in Fig. 

El 

These monotonic velocity regions are then used for counting modes. The EGD 
points are calculated by the interpolation of the E — i^n data points using the 
Vgrp- The details for calculating the energy nodes on the EGD is shown in Fig. |4] 
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As a by-product of the calculation the carrier velocity is also obtained which can 
be used for other calculations such as the mean free path. 
5. The modes from each of the ID E — Jfy are then integrated over all the Kj_ and 
divided by the unit-cell area (Auc) to obtain the total DOM. 

The present algorithm has the advantage that the original E — K^^ can be sparse 
compared to the energy grid on which the DOM is calculated since the monotonic 
regions of i5 — A'y allow to interpolate the dispersion data to be used for DOM and 
velocity calculations. The strength of the algorithm to obtain the DOM for different 
K-grid densities is shown in Fig. |5] Since the algorithm can work on sparse dispersion 
data, the time to obtain the total modes is also reduced. 

The present algorithm assumes that all the K vectors for a given energy dispersion 
are orthogonal. This assumption has both advantages and disadvantages. Since the K^^ 
and Kj_ are separated, this allows for parallel computation of modes for each K±^ set. 
This leads to computational speed-up. This aspect of the algorithm can be inspected as 
a future work. The interpolation in the E(k) is always done along ATy but not along K^. 
This puts a limitation on the sparsity of the Kj^ grid. A very sparse grid will result 
in erroneous DOM calculation. The over-all E(k) should not be too sparse either such 
that the original features of the dispersion are lost. In that case velocity interpolation 
will give erroneous results. The sensitivity of the TE results on the K-grid sparsity and 
the compute time for DOM are discussed in Sec. |3.1| and Sec |3.2[ respectively. 

Apart from the general steps adopted for the calculation of the DOM for both elec- 
trons and phonons, some special care in selecting the energy ranges for both, involved 
in Eq. (|6| and ([sjl, must be taken. 

2.3 Energy range: electron transport 

Real materials are characterized by many different electronic bands. However, not all 
these bands contribute to the electron transport and an energy range around the Fermi- 
level (Ep) needs to be selected carefully. To obtain an expedient but good approximate 
solution the energy cut-offs {Emin, Emax) are chosen such that the integral values for 
the transport parameters ( Eq. (|3|-([6|) do not show any variation. The bounds for the 
energy grid (Eq. |6| of the DOM is obtained as follows, 

Emax ^ Ec or min[max[E{K\\)\fKj^\] (10) 
Emin ~ Ev or max[min[E{K^^)'^Kj_]], (11) 

where min{max) represent the minimum (maximum) value in a given numerical 
array. The terms Ec and Ev define the conduction band minima (CBM) and the valence 
band maxima (VBM), respectively as shown in Fig. |6] Our calculations show that in 
order to obtain correct results, the Ep value can vary between the following limit, 

Emin + WkBT <Ep< Emax - WkeT (12) 

where T is the temperature. The choice of lOfc^T is chosen since this gives a good 
range where the integrals involved in the TE parameter calculations become invariant 
to the choice of energy grid as shown in Fig. [7] 
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2.4 Energy range: lattice transport 

The lattice kernel calculations do not depend on any kind of Fermi-level. Unlike the 
electron bands, the phonon bands are always within a fixed energy range with a vary- 
ing number of sub-bands for different dimensional structures |311I32| . Also there is no 
negative phonon energy dispersion in stable semi-conductor structures [33], hence the 
energy grid of the DOM always contains positive values. The energy limit for the lattice 
kernel is chosen as follows, 

^min = or User defined (f3) 
Omax or User defined, (f4) 

where Qmax is the maximum energy limit of the phonon dispersion. 



2.5 Transmission calculation 

For ballistic transport of electrons or phonons the transmission (T(E) ) of all the modes 
is 1. However, in reality carriers undergo a lot of scattering which depends on the 
dimensionality of the system, doping, temperature, etc. This reduces the transmission 
of the modes below 1. For a conductor of length L, T{E) is given by |17II23) . 

where < \[E) > is the carrier mean free path (MFP) obtained by the summing 



over all allowed k points at energy E. When L » MFP (diffusive limit) then Eq.(f5l 
can be approximated as, 

UE) . (f6) 

All the scattering mechanisms present in a system are lumped in the 'mean free 
path'. The energy dependence of the MFP can be broadly classified into two categories, 
(i) constant MFP (no energy dependence), and (ii) energy dependent MFP. For some 
common scattering mechanisms like ionized impurity, acoustic phonon, etc, < \{E) > 
can be expressed in a power law form as < \[E) >= \q[E / {kgT)Y , where E is the 
kinetic energy of the carrier, Aq is a constant and 'r' is a characteristic exponent 
describing a specific scattering process [TT]. 

In most of the BTE calculations the scattering time {jscat) is used instead of the 
MFP. Again for Tscat the energy dependence are of two types, (i) energy independent 
(constant Tscat), and (ii) energy dependent. The constant Tscat case is physically hard 
to justify since it is well known that particles scatter to/from different energy states at 
a different rate [T7] . The connection of the scattering time to the MFP is given as [T7] , 

Y,K^\\iK,E)-Tscat{K,E) 

<X{E)>=2- ' (17) 

Here the summation is over all the K states at a given energy E. If the scattering 
time is assumed isotropic in K then the MFP is given as, 
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In the present algorithm, the MFP can be calculated efficiently since the velocities 
are already obtained as a by-product during the DOM calculation. In the case of a 
constant scattering time, the energy dependence of the MFP is derived purely from 
the electronic or phonon energy dispersion. 



3 Results and Discussion 

In this section we provide the results on the dependence of calculated TE parameters 
on the sparsity of the K-grid using the algorithm. The timing analysis is provided to 
give an idea about the total speed up with K-grid reduction and which part of the 
calculation consumes the maximum compute time. Also the comparison of the DOM 
and TE parameters, calculated from the algorithm, with analytical expressions and 
realistic material dispersions are provided. 



3.1 Sensitivity analysis: How robust is the algorithm? 

To understand the strengths and limitations of the algorithm, K-grid sensitivity tests 
are performed using parabolic E(k) dispersions for 3D, 2D and ID cases. The param- 
eters used for the generation of the parabolic bands are shown in Table [2] 

The TE parameters like G and S are calculated using the parabolic E(k). The 
k-grid density variation introduces numerical error in the S and G calculation. The 
percentage error in the maximum power factor (PFmax) is related to the error in S 
and G by the following relation, 

APF^nax _ , AS AG 

where AB/G are the variations in Seebeck coefficient and electronic conductivity, 
respectively. The final fluctuation in the PFmax depends on the sign of AS/G. However, 
all the fluctuation plots are shown for the absolute value of the errors. 

To start the sensitivity analysis, first a base K-grid is chosen. A K-grid with 100 
points in each direction (-Tr/ag to n/ao) is found to be sufficient to provide stable 
results. A finer K-grid does not change the final results by more than 0.5% for any of 
the calculated TE values. Three different types of studies were performed to determine 
the sensitivity of the algorithm to reduction in K^^ and Kj_, 

— Case A: Keep fixed at 100 grid points and reduce K^^ down to ~ 60%-80%. 

— Case B: Keep K^^ fixed at 100 grid points and reduce Kj^ down to ~ 60%-80%. 

— Case C: Reduce both K^^ and K± down to ~ 60%-80%. 

For the sake of brevity only the 2D TE error analysis results are shown. Other 
dimensions (3D and ID) show similar results and the outcomes are similar. For the ID 
system only case C is applicable since in these systems only ify is the free momentum 
direction. The other two directions are geometrically confined as shown in Table [l] 



8 



Table 2 Parameters used for the generation of parabolic bands 



Structure 






■m^2 


Ec 


En 


ao 


(Dim) 


xmo 


xmo 


xmo 


eV 


eV 


nm 


Wires (ID) 


1 






0.2 


-0.2 


1 


Films (2D) 


1 


1 




0.2 


-0.2 


1 


Bulk (3D) 


1 


1 


1 


0.2 


-0.2 


1 



Sensitivity Analysis: Case A 

The reduction in J-f|| down to 60% results in less than 1% variation in S and ~6% 
variation in G as shown in Fig. [8|i. The corresponding fluctuation in the PF is around 
5% as shown in Fig. [sja. The Fermi-level, at which the maximum PF is extracted, 
however remains unchanged. The fluctuation in the TE parameters arises only from 
the fluctuation in the DOM. Thus, the present DOM calculation method is quite robust 
to reductions in ify given the K^^ has good mesh density. 

Sensitivity Analysis: Case B 

The reduction in K^^ down to 60% results in less than 2% variation in S and ~12% 
variation in G as shown in Fig. |9]j,. The maximum fluctuation in PF is around 10% as 
shown in Fig.|9|3. The Fermi-level {Ep) at which the maximum PF is extracted shows a 
maximum variation of ~2.5%. In this case, the fluctuation in the TE parameters arises 
from the fluctuation in, (i) the DOM, and (ii) the Ep. The present DOM algorithm is 
sensitive to variations in Kj^_. 

Sensitivity Analysis: Case C 

The reduction in all the K points down to 60% results in less than 5% variation in S 
and ~13% variation in G as shown in Fig. |10^ . The maximum fluctuation in PF is 
around 10% as shown in Fig. |10| 3. The Ep at which the maximum PF is extracted 
shows a maximum variation of ~2.5%. Thus, the fluctuation in the TE values arises 
from the fluctuation in, (i) the DOM, and (ii) the Ep. This case has almost similar 
K-grid sensitivity as case B, again showing that the present DOM algorithm is sensitive 
to variations in K^^. 

3.2 Timing analysis 

The present algorithm can calculate the TE parameters within reasonable error limits 
as shown in the previous section. Another obvious question that arises is how much 
computational speed-up can be achieved. The time to calculate the DOM for the three 
cases presented in the previous section is analyzed for 3D, 2D and ID structures. 

As the K density along all the directions is reduced, the speed up for each dimension 
is different. For the 3D system, the time required goes up with total number of K-points 
{NK) with a power of 1.46 [NK'^'^^). For the 2D case the power law is NK^ '^^ and for 
ID case the time taken is almost constant (in the given NK range). The algorithm takes 
roughly 900 seconds for 1 million K-points ( 100 x 100 x 100) for 3D case on nanoHUB.org 
workspace |34]. For the 2D case the time taken for ten thousand K-points (lOOx 100) is 
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Table 3 Summary of sensitivity and timing analysis 



Dimension 


K-grid 
direction 


Max. K 
reduction (%) 


S 


Max. 
G 


3rror (%) 
PF 


Ep 


Good E{K) 
sparsity 


DOM speed up 
(60% K red.) 


3D 


^11 


80 


-4.41 


3.02 


-5.65 


4.14 


<70% 


1.02X 


3D 




64 


-5.3 


-4.45 


-11.49 


1.02 


<50% 


8.4x 


3D 


All-K 


66 


-5.3 


15.2 


5.9 


5.2 


<50% 


3.6x 


2D 


^11 


60 


-0.4 


6 


5 





<60% 


l.lx 


2D 




60 


-0.32 


-9.30 


-9.95 


1.18 


<50% 


1.8x 


2D 


All-K 


58 


4.2 


-13.3 


-4.4 


1.18 


<50% 


1.5x 


ID 


All-K 


80 


4 


13.2 


20 


3.8 


<70% 


1.05X 



nearly 2 seconds and for ID case the time taken is roughly 0.1 second for 100 K-points. 
All these results are shown in Fig. 

For the cases of 2D and 3D, the algorithm requires different compute times with 
reductions in K-points along both the transport and the transverse direction. The 
compute time for the DOM {tooKl) is roughly independent of the K-point reduction 
in transport direction for both the 2D and 3D cases (Fig. 12 1. However, for a 60% 



reduction in K_[^, the 2D case shows a tuOM speed up factor of ~2 (Fig. 12 a). While 
for the 3D case, a speed up factor of 6 is observed (Fig. [12] b). A reduction in all 
K-points along all directions show a similar speed up (Fig. |12|. Thus, the present 
algorithm shows a good speed up with point reduction. 



3.3 Discussion: Algorithm aspects 

The TE sensitivity analysis and tjjoM speed up reveal that the algorithm to calculate 
the DOM is more sensitive to the Kj_ points compared to the K^^ points. A reasonable 
reduction in Kj^ must be chosen in order to optimize the compute time and to obtain 
reasonably stable TE parameter values. A summary of all the analysis is provided 
in Table |3] This table also provides the limits for reduction in K points in the E(k) 
data-set to obtain TE parameters within a 10% error margin. In most of the cases a 
50% reduction in K-points is easily achievable without a big penalty on the calculated 
TE values. The sensitivity analysis presented here is for parabolic bands, however, 
the general features of the algorithm remain quite similar even for the dispersions of 
real materials which are more complex than parabolic bands. Similar conclusions are 
obtained for the phonon dispersions. 



3.4 Calculation of the TE parameters 

The final verification of the algorithm is done by calculating the TE parameters for (i) 
the parabolic bands in 3D, 2D, and ID cases, and (ii) bulk Bi2Te2,. 

Parabolic Bands 

The parameters used for the electronic energy dispersion are shown in Table |2] For all 
dimensions the number of energy points in the DOM (EGD) is kept constant at 500. 
The analytical results for the modes and TE parameters are obtained from Ref. [17] 
and [18]. The number of modes for all three dimensions compare very well (< 4% error) 
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with the analytical results as shown in Fig. |13| Using the modes, the TE values are 
calculated. Only the 2D case is shown in Fig. |14[ The agreement with the analjrtical 
calculations is very good with around 1% error in the numerical values. The 3D and 
ID cases also show very good agreement with the analytical calculations. Thus, the 
algorithm provides accurate results for the electronic TE parameters. 



Realistic Bands: Bulk Bi2Te^ 



As a next step of verification, the algorithm is tested for the calculation of TE pa- 
rameters for Bi2Te^. The same Tight-binding electronic dispersion [35] is used for 
the calculation of the DOM as used in Ref. [13 ■ The agreement with the published 



DOM results is very good (within 1% error) as shown in Fig. 15 1. Using the DOM, 
the S and G axe also calculated which are used to obtain the PF. The calculated PF 
again shows a very good agreement with the published theoretical result [17] as well 
as with the experimental data [36^ as shown in Fig. 15 a. The same calculation is also 



performed for the lattice thermal conductivity of bulk Bi2Te^. The phonon dispersion 
is obtained using GULP [37] as provided in Ref. [53]. The agreement of the calculated 



phonon modes with the published theoretical result [23] is very good (Fig. 16 1). Also 



the lattice thermal conductivity calculated using the method provided in Ref. |23) gives 
a very good agreement with the experimental value [3S]. Thus, the present algorithm 
provides accurate TE values for real materials too. 



4 Summary and Outlook 

An efficient algorithm to calculate the electron and phonon modes in any dimension 
is presented. The algorithm provides an efficient implementation of a TE parameters 
calculation scheme based on the Landauer's approach and will be extremely useful in 
readily and accurately evaluating the existing as well as new thermoelectric materials. 
The algorithm is sensitive to the transverse K point density in the E(k) relation both in 
terms of the final TE calculations as well as the compute time. A proper optimization 
of the K-point reduction is provided to allow for fast and accurate TE parameter 
calculations. The results from the algorithm are also bench-marked with analytical 
and real material TE parameter values. This algorithm will be useful for developing 
computer programs to evaluate the TE performance of new and artificial materials in 
the future. 
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Fig. 1 Modeling hierarchy for the thermoelectric analysis. The present work focuses on the 
calculation of the TE parameters from the energy dispersion relations as shown by the encircled 
part. 




Fig. 2 The band-counting method for calculating the DOM. (a) The lowest sub-band of the 
electronic E(k) of a 20nm X 20nm [100] SiNW. (b) The corresponding propagating modes 
M(E) associated with this single band. 
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Fig. 3 Velocity at eacli point positive half of first conduction band sub-band for a 2nm X 
2nm [100] Si nanowire (SiNW). As mentioned above, the important points to note in the given 
band are the points where sign of velocity changes. These points are indicated as PI, P2, and 
P3 and corresponding regions of interest are marked as Rl, R2, and R3 in Figure 2.3. Each 
point provides monotonic velocity range (increasing or decreasing) and calculations for DOM 
are performed separately for these ranges within a sub-band. 



E-KGrid DOME- Grid 




Fig. 4 The E — K points on the provided energy dispersion shown by dots. The energy 
grid for DOM (EGD) points are shown using crosses. These points are of two types, (i) the 
matching point to the E — K grid, and (ii) the one which requires interpolation of the provided 
E — K relationship. This interpolation (either linear or quadratic) is done in the appropriate 
monotonic energy region like iJ2,etc shown in Fig.[3] In this way the DOM is created for 
the EGD. 
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Fig. 5 (a) E{K) relation with different number of k points. Case A with 61 K points and case 
B with 13 K points. Also the two monotonic E — K regions are shown (Rl and R2) along with 
the turn around point PI. DOM calculated for the two E — K grids using (b) 50 energy points 
and (c) 12 energy points. The DOM matches 100% for all the 4 cases showing the robustness 
of the DOM calculation method. As long as the sparse E — K captures the important turn 
around points (like PI) correctly the DOM calculation algorithm obtains the correct number 
of modes. 




Fig. 6 Schematic showing the range of energy limit and the range of Fermi-level used for the 
calculation of the integral in Eq.|6]for electrons. 
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Fig. 7 Variation in the values of (a) /q E'nd (b) Ii (Eq.[6| for different values of the Fermi- level 
(Ep). When the Ep is at least lOfcsT below the Emax then the integral values show less than 
1% variation. Similar result is also obtained for the integral l2- 




Fig. 8 Impact of K\i^ point reduction on (a) S (left) and G (right) and (b) PF (left) and Ep 
(right) in a 2D structure. All the values are extracted at the maximum PF point. The has 
100 grid points. Even for 60% reduction in i^y points none of the TE values show more than 
6% variation. 
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Fig. 9 Impact of K±^ point reduction on (a) S (left) and G (right) and (b) PF (left) and 
Ep (right) in a 2D structure. All the values are extracted at the maximum PF point. The K\i^ 
has 100 grid points. For 60% reduction in K-points G shows a maximum variation of 12% and 
PFmax has variation around 10%. 
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Fig. 10 Impact of reduction in all the K-points, on (a) S (left) and G (right) and (b) PF 
(left) and EF (right) in a 2D structure. All the values are extracted at the maximum PF point. 
G shows a larger fluctuation (>10%) compared to S fluctuation (<4%) which also reflects in 
the PFmax fluctuation. 
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Fig. 11 DOM compute time for ID, 2D and 3D parabolic bands. The number of K points 
are reduced along all the K-directions equally. The 3D case takes the maximum time due to 
higher number of K-points, followed by the 2D and the ID case. 
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Fig. 12 DOM compute time (toOM) for the three types of K-point reduction for the (a) 2D 
structure and (b) 3D structure. For both the cases the compute time is almost independent of 
/fli reduction (brown Une). reduces with K± point reduction. For each point 5 compute times 
are averaged. 
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Fig. 13 Comparison of the numerical modes calculation using the algorithm with analytical 
modes calculation using parabolic bands with m* = mo (from Ref \18\ ) for (a) 3D, (b) 2D and 
(c) ID structure. The steps in the 2D case appear due to the sparse energy grid chosen. 
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Fig. 14 Comparison of the numerical calculation with analytical expression for effective mass 
from Ref 18j for a 2D system for (a) Conductance (b) Seebeck Coefficient and (c) Thermoelec- 
tric Power Factor at 300K. The numerical results compare within 1% to the analytical values. 
The parameters used for the parabolic bands are provided in Table [2] 
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Numerical 




Fig. 15 Comparison of (a) DOM calculation and (b) Power factor at 300K for Bi2Te3, using 
the algorithm, with the theoretical calculations reported in Ref. 17 and experimental results 
from [36[ . The electronic energy dispersion for bulk Bi2Te-i is obtained using the sp^d^s* 
tight-binding model [35) . The PF matching for Bi2Te3 is obtained assuming a constant mean- 
free-path of 18, 4 nm for conduction and valence bands respectively as reported in Ref. |17| . 
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Fig. 16 (a) Comparison of the bulk Bi2Te^ phonon modes calculated using the algorithm and 
theoretical value reported in Ref. |23) . The bulk phonon dispersion is obtained using GULP 
|37| . (b) Comparison of the simulated and experimental 1381 thermal conductivity for Bulk 
Bi2Te3 from 50 to 500K. The phonon scattering mechanisms considered here are outlined in 
detail in Ref. [23]. 



